#include <complex.h>
#include <math.h>
static double _Complex poly_val2(double p[], int order, double omega)
{
    int i;
    double _Complex z, sum = 0;

    sum = 0.0;
    z = cos(omega) - I * sin(omega);
    for (i = order; i >= 0; i--)
    {
        sum = sum * z + p[i];
    }
    return sum;
}

/**
 * 计算数字滤波器的频率响应
 * num 是数字滤波器的分子多项式系数
 * den 是数字滤波器的分母多项式系数
 * num_order 是分子多项式的阶数
 * den_order 是分母多项式的阶数
 * sign = 0 时，x_out 为频率响应的实部， y_out 为频率响应的虚部
 * sign = 1 时，x_out 为频率响应的模， y_out 为频率响应的幅角
 * sign = 2 时，x_out 为以 dB 为单位的频率响应， y_out 为频率响应的幅角
 * len 为，频率响应的取样点数
 */
void gain(double num[], double den[], int num_order, int den_order,
            double x_out[], double y_out[], int len, int sign)
{
    int k;
    double omega;
    double _Complex h;
    for (k = 0; k < len; k++)
    {
        omega = M_PI * k / (len-1);
        h = poly_val(num, num_order, omega) / poly_val(den, den_order, omega);
        switch(sign)
        {
        case 1:
            y_out[k] = carg(h);
            x_out[k] = cabs(h);
            break;
        case 2:
            y_out[k] = carg(h);
            x_out[k] = 20 * log10(cabs(h));
            break;
        default:
            x_out[k] = creal(h);
            y_out[k] = cimag(h);
            break;
        }
    }
}
